First-Principles Based Matrix-Green's Function Approach to Molecular Electronic 

Devices: General Formalism 

Yongqiang Xue * 

Department of Chemistry and Materials Research Center, Northwestern University, Evanston, IL 60208 and School of 
Electrical and Computer Engineering, Purdue University, West Lafayette, IN 47907 

Supriyo Datta 

School of Electrical and Computer Engineering, Purdue University, West Lafayette, IN 47907 

Mark A. Ratner 

Department of Chemistry and Materials Research Center, Northwestern University, Evanston, IL 60208 

(February 1, 2008) 

Transport in molecular electronic devices is different from that in semiconductor mesoscopic 
devices in two important aspects: (1) the effect of the electronic structure and (2) the effect of the 
interface to the external contact. A rigorous treatment of molecular electronic devices will require 
the inclusion of these effects in the context of an open system exchanging particle and energy with 
the external environment. This calls for combining the theory of quantum transport with the theory 
of electronic structure starting from the first-principles. We present a rigorous yet tractable matrix 
Green's function approach for studying transport in molecular electronic devices, based on the Non- 
Equilibrium Green's Function Formalism of quantum transport and the density-functional theory 
of electronic structure using local orbital basis sets. By separating the device rigorously into the 
molecular region and the contact region, we can take full advantage of the natural spatial locality 
associated with the metallic screening in the electrodes and focus on the physical processes in the 
finite molecular region. This not only opens up the possibility of using the existing well-established 
technique of molecular electronic structure theory in transport calculations with little change, but 
also allows us to use the language of qualitative molecular orbital theory to interpret and rationalize 
the results of the computation. We emphasize the importance of the self-consistent charge transfer 
and voltage drop on the transport characteristics and describe the self-consistent formulation for 
both device at equilibrium and device out of equilibrium. For the device at equilibrium, our method 
provides an alternative approach for solving the molecular chemisorption problem. For the device 
out of equilibrium, we show that the calculation of elastic current transport through molecules, both 
conceptually and computationally, is no more difficult than solving the chemisorption problem. 

I. INTRODUCTION 

There has been significant progress in exploring the concept of molecular electronics in recent years, due to the 
advancement of techniques for characterizing and manipulating individual moleculesETEl. The fact that useful .devices 
can be built on the basis of individual molecules, as demonstrated recently by several research groups l3~1i3, has 
generated wide-spread interest in this new technology. In order to fulfil the true promise of molecular electronics, it 
is essential to have a thorough understanding of the electronic and transport processes at the single molecule level. 
This paper represents an attempt to put our understanding of electronic transport through individual molecules on a 
firm theoretical basis starting from the first-principles. 

Traditionally electron transport phenomena are studied in the context of bulk semiconductor devices, the theoretical 
description of which is largely built upon two premises: (1) the effective-mass equation and (2) the Boltzmann 
Transport Equation (BTE)0J. The effective-mass equation subsumes the effect of the background periodic lattice 
potential into an effective Hamiltonian so that the electrons can be considered as particles of effective mass m* 
in some applied field. The semi-classical nature of the electronic motion in the bulk devices, on the other hand, 
allows us to describe the distribution of electrons in response to the applied fields and various scattering mechanisms 
through the solution of the BTE equation, in much the same way as that of classical point-like particles. Within this 
approach, the quantum-mechanical effect only comes in through the calculation of band structures and the various 
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carrier scattering rates, which provides the input to the solution of the BTE or its approximate versions such as the 
drift-diffusion equation. As a result, the study of transport phenomena is effectively decoupled from the study of the 
electronic structures. 

The investigation of quantum mechanical transport has begun to flourish during the past two decades, largely 
due to the advancement in lithographic techniques ..which has allowed routine fabrication of submicron features in 
artificially tailored semiconductor heterostructuresHEd. Two quantum mechanical effects distinguish such mesoscopic 
devices from bulk devices, reflecting the wave-particle duality of electron. Onc_is the quantization of electronic 
charge which evidences itself in coulomb blockade and single-electron transistorsEj. The other is the preservation 
of quantum phase coherence over a length with size comparable to one of the device dimensions and the resulting 
energy quantization of confined electrons, which leads to the observation of conductance quantization in transport 
through a narrow constriction (quantum point contact) and negative differential resistance in double-barrier tunneling 
structures (resonant-tunneling diode). Transport in these quantum semiconductor structures is now determined by 
the scattering property and occupation of the electronic eigenstates under an appropriate external potential. On the 
other hand, electrons in such quantum semiconductor structures are usually confined in one or two directions. The 
confinement potential, often the electrostatic potential due to a nearby gate or the band discontinuity across the 
heterostructure interface, varies slowly on the atomic scale. So the electrons can still be described satisfactorily as 
free carriers using the single-band effective mass approximation and most of the import ant p henomena in mesoscopic 
transport can be understood without worrying about the complex electronic structures! 20 ^ 1 !. 

The situation changes dramatically when we consider molecular electronic devices formed by sandwiching a chemi- 
cally synthesized molecule between two large (on the molecular scale) metallic electrodes. At the molecular scale, the 
simplicity associated with the effective-mass approximation breaks down and the electronic structure of the system 
has to be taken explicitly into account. The quantum mechanical scattering problem involved is now the scatter- 
ing of electrons under the potential of the atomic nuclei and the potential due to other electrons. Understanding 
transport in such molecular devices therefore requires a knowledge of the microscopic electronic structure of the 
electrode-molecule-electrode system under both zero and finite bias voltages. 

The reduction to molecular scale also brings in another complication in transport study as compared to the meso- 
scopic systems: treatment of the interface to an external contact. In mesoscopic transport, the details of contact 
are often not important. The measuring electrodes, taken as infinite electron reservoirs, can either be simulated by 
reflectionless semi-infinite leads with simple confinement potential at the interface, or only come into the theoretical 
formulation as an appropriate boundary conditiorEj'Eil. This is no longer true when the device is of molecular di- 
mension. Since the electrodes can have atomic structures on the surface whose dimensions can be comparable to the 
molecule, the usually well-defined boundary between the active device region and the contact region is blurred. The 
interface to the external contact becomes an integral part of the device and the measured electrical characteristics will 
depend on the details of the atomic arrangement of the contact. Moreover, the electronic and structural propectiss-pf 
the molecule could be modified by the bonding to the measuring contact, bringing in additional complicationsE3~E3. 

In summary, molecular electronic devices are different from their mesoscopic counterparts in two important aspects: 
(1) the effect of the electronic structure and (2) the effect of the interface to the external contact. Since the molecule 
can freely exchange energy and electrons with the electrodes, a rigorous treatment of molecular electronic device 
can only be achieved including these effects in the context of an open system. As a result, a successful modeling 
of molecular electronic devices in general calls for combining the theory of quantum transport and the theory of 
electronic structure starting from first-principles. 

There have been numerous theoretical works on transport through individual molecules describing the electronic 
structure at different levels. Early works have focused on understanding the fundamental mechanis m s o f-transport 
in the molecular-scale and developing simple theory for the explanation of experimental results EjoEH^. These 
works thus were centered on model Hamiltonian or semi-empirical theories (notably the Extended-Hiickel-Theory 
of organic moleculesE-3 and 7r-electron tight-binding theory of carbon nanotubesE3) . Due to the interdisciplinary 
nature of molecular electronics, different methods have been used reflecting-the. authors' own background which are 
essent ia l ly -all-equivalent to the Landauer formula of mesoscopic transportE3~E3 as used extensively in our previous 
work«§M. Such works have provided useful insights into the factors governing transport through individual 
molecules. However, the usefulness of this approach is limited by its incapability of providing an accurate description 
of the electronic structure of the molecule and the metal surface involved. Even if a good parameterization exists 
for the molecule and the bulk separately, the charge transfer and the resulting self-consistent charge and potential 
relaxation upon the formation of surface and adsorption of the molecule is difficult to treat unless drastic assumptions 
are made or additional parameters are introduced. In addition, it is difficult to take into account coupling between 
electrons ancLmplecular vibrational modes which is expected to play an increasingly important role as the molecular 
size increasesEjEJ. I-, .— . 

More recently, Lang and coworkersc3~c3 have presented self-consistent studies of both the conductance and current- 
voltage characteristics of atomic and molecular wires using the jellium model of a metal surface and the local- 
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density- approximation (LDA) of density functional theory (DFT) with plane- wave basis set. The calculation proceeds 
by recasting the Kohn-Sham equation of the electrode-molecule-electrode system into a scattering riorm using the 
Lippmann-Schwinger equation and solving the wavefunctions of the scattering state self-consistentlycHl. The current 
is obtained by summing over the contribution of the scattering states which are occupied according to the Fermi 
distribution in each electrode, following the spirii—of Landauer theory. Similar approaches have also been used by 
Tsukada and coworkerscH and Guo and coworkers^. 

The jellium model has an appealing and physically reasonable simplicity: all the complexities of the bulk band 
structure are simply ignored, the effect of the substrate persists only in providing a continuous energy spectrum and 
the only inhomogeneity left in the metal region is the essential one-dimensional inhomogeneity of the surface itself. 
However, the jellium model is known to be deficient in describing the electronic density of states and charge densi ty , in 
the region perturbed by the absorbed molecule even for sp-bonded simple metals for which it is more successfult^E3, 
and it cannot answer questions regarding the effect of adsorption geometry and surface relaxation. In addition, it is not 
applicable to semiconductors and gives serious quantitative error when applied to noble and transition metals where 
the bonding at the surface is more directional. The use of the plane-wave basis set and the wavefunction formulation 
also make it difficult to treat larger molecular systems, to improve the description of the electronic structure and to 
include other scattering mechanisms. 

The objective of this paper is to present a rigorous yet tractable self-consistent matrix Green's function (MGF) 
approach for studying transport in molecular electronic devices, through which many limitations of the previous 
works can be overcome rather conveniently using existing computational techniques. In particular, this allows us to 
use existing quantum chemical methods for transport calculation with minimum change. The method is based on the 
Non-equilibrium Green's Function Formalism of quantum transport^] and the density-functional-theory of electronic 
structureEJ. We also use the technique of expansion in a finite set of local orbital functions. 

The Non-Equilibrium Green's Function Formalism has proved to be a powerful and formally rigorous approach to 
studying quantum transport phenomena in nanodevices, where the devices can be characterized as a small active 
device region connected to electron reservoirs via non- interacting leads. Arbitrary interactions in the device can be 
includeoH'E-X For nop-iisteracting electron systems in the coherent limit, it reduces to the familiar Landauer theory of 
mesoscopic transport! 5 tPi. The formalism is fairly complete, and it has been applied successfully to studying quantum- 
mechanical transpGaJji^jhenomena in semiconductor mesoscopic systems such as resonant-tunneling diode and single- 
electron transistorsEjOHjjjsually in combination with a model Hamiltonian and more recently using semi-empirical 
tight-binding formulationso. For molecular electronic devices, within the Born-Oppenheimer approximation, the 
problem is to calculate the Green's function of the interacting electron system under the potential of the given 
configuration of atomic nuclei and external fields. 

The study of molecular electronic devices is greatly facilitated by recognizing the fact that due to the metallic 
screening in the electrodes, the charge and potential perturbation induced by the adsorption of the molecule extends 
over only a finite region into the electrodes, practically the region enclosing the surface metal atoms closest to the 
molecule. The charge and potential distributions beyond this region are the same as that of the bimetallic interfaces 
without the molecules. By expanding the wavefunction in terms of the finite set of local atomic orbital functions, 
the matrix Green's function approach allows us to take full advantage of this spatial locality by separating the device 
into the "extended-molecule" (which includes the molecule itself and the surface atoms perturbed by the molecule) 
and the electrode region, the description of each can be treated and systematically improved independent of each 
other. The effect of the external contact can be included rigorously as a self-energy operator, which depends only on 
the charge and potential distribution outside the "extended molecule" and needs to be computed only once, thereby 
the computation can be focused on the physical processes happening in the finite "extended molecule" region. This 
not only opens up the possibility of using the full repertoire of molecular electronic structure calculation for studying 
transport in molecular electronic devices, but also allows us to interpret the result of computation using the simple 
picture of bonding and orbital interactions in molecules, which is necessary for the computation to be useful in terms 
of understanding rather than pure numbers. Although we will focus on the density-functional theory of the interacting 
electron system, the formalism doesn't depend on the particular description of the electronic structure. 

In .matrix form, the present formulation is equivalent to previous works using semi-empirical tight-binding formula- 
tionsEaO. The tight-binding formulation provides a conceptually and computationally simple framework for studying 
transport in systems where a quantum mechanical description of the electronic structure is necessary but the charac- 
teristic length scale is much larger than inter-atomic spacing so atomic-scale details are not needed. Example systems 
include layered semiconductor devices and long carbon nanotubes. The present formulation keeps the simplicity of 
the tight-binding approach while putting it on a firm theoretical basis. The explicit use of the local basis function in 
the real space also allows much more detailed understanding of the physical processes through quantities such as the 
spatial distribution of chargej-arui current densities, which play the fundamental role in density-functional theory and 
its time-dependent extensionElfLJ. In addition, the ambiguities associated with the semi-empirical formulation, such 
as the orthogonality of the basis functions used in such formulation and the treatment of self-consistent charge-transfer 
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effected can be avoided. When the potentials and the integrals are calculated approximately rather than from first 
principles, it ca» also serve as the basis of improved self-consistent tight-binding approaches capable of handling much 
larger systemsEl 

The remainder of this paper is organized as follows: We discuss briefly the use of Non-Equilibrium Green's Function 
and density-functional theory for modeling transport in molecular devices in Sec. II. The matrix Green's function 
method is described in Sec. III. The self-consistent formulation for device at equilibrium is given in Sec. IV, where 
we show that our method provides an alternative and generalizing approach to the familiar chemisorption problem in 
surface physics. The self-consistent formulation for a device out of equilibrium is given in Sec. V. Finally we devote 
Sec. VI to conclusions. We use atomic-units throughout this paper unless otherwise noted. 



II. NON-EQUILIBRIUM GREEN'S FUNCTION APPROACH FOR MODELING MOLECULAR 

ELECTRONIC DEVICES 



A. Non-Equilibrium Green's Function Formalism 

As current flows, the device is driven out of equilibrium. For systems out of equilibrium, the Green's function 
approach can be-developed following the same procedure as the equilibrium case, by defining the contour-ordered 
Green's functiorSl^ G(l,l') = -i(T c [ip H (l)^ + (l')]) (for details, see Haug and JauhcH). The Green's functions 
involved, besides the retarded and advanced Green's function, 

G R (i, 1') = -iQ(h - tiOttV'* (i), V£(i')}>, 

G A (i,i')=ie(t 1 ,-*i)({^(i),vJ(i')}), (2.i) 

include the correlation function (or the "lesser" Green's function) which is the central quantity in this formalism 

G<(f,i;rV') = +i^ + (r,t)ip(r',t')) (2.2) 

Any physical observable can be obtained from G < (r, t; f", t') and its transformations. For example, the quantities of 
most interest to us, the charge density n(r,t) and the current density j(r,t), are determined as following: 

n(r,t) = ^+{r,t)^{r,t)) = -iG< (r, t; r, t) (2.3) 

and 

j(r,t) = 1 Jhn (V - V)(^+(f,t)^,t')) = I lim_(V' - V)G<(r,t;rV) (2.4) 



For steady state, which we consider here, the Green's functions depend only on the time difference t — t , which we 
can Fourier transform to energy. The resulting equations of motion (EOM) of the Non-Equilibrium Green's function 
are the Keldysh-Kadanoff-Baym equatior£3 EJ: 

{E - [~V 2 + V ext {r)]}G R {r,7;E) - J dr"Y, R (r, r"; E)G R (r" , r 1 ; E) = S(f- ?), (2.5) 

and 

{E- [-\^ 2 + V ext {r)]}G<{r^ ] E) - J dr"S«(r,r" ; £)G<(r" ,f ; E) 

= J d^^ < {r^'- J E)G A {^'^-E) (2.6) 

or its equivalent integral formulation: 

G< (r, r'; E) = J dr" J dr*"G R (r, r"; £)£<(^', f"; E)G A {?" , r<; E) (2.7) 

where we assume that the non-equilibrium term in the Hamiltonian is incorporated into a one-body external potential 
V ex t- The interactions are contained in the self-energy operators S[G] which can be obtained-, systematically from 
many-body perturbation theory following the same procedure as their equilibrium counterpart!^!. 
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By solving the KKB equation for G < , we will be able to calculate the current distribution within the device given 
the Hamiltonian of the system. However, in general, they involve the evaluation of the full correlation function and 
the retarded Green's function in the presence of tunneling into the leads. The usefulness of the above expressions 
then depends on whether we can devise practical calculation schemes for the single-particle Green's functions. For 
molecular-scale devices, the problem is to determine the state of the interacting electron system under the potential of 
the given atomic nuclei. By separating the term describing the electron interactions into the classical Coulomb part 
and the exchange-correlation part, we can write down the Keldysh-Kadanoff-Baym equation in the following form: 

{E-[-\v 2 + V ext (r) + J d7^^}}G R (ry;E) 
dr'z R (r,r';E)G R (r',r;E) = S(r-r), 
G < {r,r;E) = I dr" I dr" G R (r, r'\ E)T, < [f\ r"; E)G A (r" ,r ;E) (2.8) 



where the external potential V ex t is the summation of the ionic potential and the potential that drives the system out 
of equilibrium (see the discussion in sections IV and M later) . 



B. Density- Functional Theory as the Compromise 

For transport through molecular scale devices, the system is characterized by the highly non-equilibrium distribution 
in the device region, which is driven by contact to two large reserviors with different electrochemical potential. A 
truly first-principles treatment of the electronic processes in such non-equilibrium system can onliz_.be based on the 
non-equilibrium version of the many-body, perturbation theory such as the quasi-particle theory_3 or perhaps the 
time-dependent density- functional, theory--. Such a first-principles theory doesn't exist yet. Instead, we will work 
with density-functional theory____l. 

The theory of DFT is based on the study of systems in their ground state or in thermodynamic equilibrium 
corresponding to the grand canonical ensemble, where the variational formulation of quantum mechanics can be 
used. For transport through non- interacting systems, the problem can be solved by working in the scattering state 
representation, since electrons coming from different electrodes are in separate thermodynamic equilibrium at different 
electrochemical potentials. This is the essence of the Landauer theory of quantum transport. In this case, it may 
be justified to use density-functional theory to calculate the scattering state wavefunctions which correspond to a 
two-component system, each characterized by a well-defined thermodynamic ensemble. However, we are not aware 
of any formal .pxpoLof this. In fact, this has been the basis for the use of DFT to calculate current and conductance 
in the pastciFll_]'_Zr_a. For the interacting electron system which DFT attempts to describe, the nonlinear current- 
voltage characteristic can be obtained only by doing perturbation theory in the part of the Hamiltonian that drives 
the system out of equilibrium. For these reasons, the use of the DFT formalism, with probably the exception of 
its time-dependent extension!-., in transport calculations can only be taken as qualitative in principle, although this 
doesn't preclude its quantitative success in practice. 

From here on, we will work with density functional theory in our modeling of molecular devices, with the un- 
derstanding that we approximate the true self-energy operator for exchange-correlation interaction with the DFT 
description of exchange-correlation potentialE-l. We will treat the usage of DFT in a more qualitative sense in that 
it has a well-defined physical basis, from which we know where it can be expected to succeed and where it may fail 
and why. Therefore we will put more emphasis on using the results of such calculation for qualitative understanding 
whenever possible. These results may serve as the basis for further improvement as our understanding of molecular 
devices progresses. The reason for this choice is therefore practical, rather than fundamental. 



III. THE METHOD OF MATRIX GREEN'S FUNCTION 



Our starting point is the Keldysh-Kadanoff-Baym equation (Eq. (2.£)). Remember that we approximate the self- 
energy operator £ by the DFT description of the exchange-correlation potential V xc (f, f") which is an energy inde- 
pendent real operator, but doesn't need to be local. In particular, we will only assume V xc to be local in the electrode 
region. It can take any non-local form in the "extended molecule". There is no "lesser" self-energy operator S < 
associated with V xc . This is in contrast with the true quasi-particle theory where £ is in genecaLnon-Hermitian and 
the corresponding £< represents the scattering rates or "life-time" of the quasi-particle stateE_H_-. 
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Given the exchange-correlation potential V xc (f, f) , we can define the single-particle wavefunction ip^r): 
-1/2V 2 + V ext (r) + f 



The retarded Green's function G R is related to ipfj,(r) through the spectral representation: 



(3.1) 



(3.2) 



where E + — lim ( s_ + o+ E + iS. We expand the wavefunction in terms of a finite set of local orbital functions: 

N 



(3.3) 



where 4>i{r) are atom-centered and decay rapidly to zero away from the corresponding atomic center. We use the 
symbol = to indicate the above expansion is exact only for a basis set that is complete. Besides the approximation 
of single-particle theory, this is the only approximation involved in our matrix Green's function method^. Note the 
choice of the basis set can be different in different parts of the system, reflecting the different nature of the electronic 
states in each of them. The Schrodinger-type equation can be transformed into a generalized eigenvalue problem: 



5> 



where Sij is the overlap matrix, 
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^r0*(r%-(r) 



and 



Hi. 



d s r<t>*{r)[~V 2 + V ext (r) + I d 3 / 



P? 



d 3 r / JPr f tf(f)V lcc (?,r)]<l> j (?) 



(3.4) 



(3.5) 



(3.6) 



Instead of solving the above equation, we substitute Eq. (3.3) into Eq. (3.2) and write down the Green's function in 
terms of basis function <bi : 



E+ - e 



T-,3 



where 



G«jiE) = Y, 



E+ - Si 



(3.7) 



(3.8) 



It is straightforward to prove that Gfj (E) satisfies the following matrix equation: 

J2(E+S lk -H tk )G«(E) = S t3 



(3.9) 



The above equation is an infinite matrix equation involving orbital functions centered around every atom of the 
molecular device. 

Since we are interested only in the Green's function in the "extended molecule" , we divide the system into three 
parts correspond to the left electrode, the right electrode and the "extended molecule" and write down the above 
equation in block matrix notation: 



/ E + Sll — Hll E + Slm — Hlm E + Slr — Hlr 



III 

E + S M l - H ML E + S M m - H MM E + S M r - H MR \ x | G ML G MM G MR | = | Imm 
y E + Srl ~ Hrl E + Srm — Hum E + Srr — Hrr 



r<R nR nR. 

^LL U LM ^LR 



'ML ^ MM ^MR 
r-iR r~<R (~*R 
^RL ^RM ^RR 



(3.10) 







Irr 
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The short-range of the local orbital basis means that for any reasonable molecule size and electrode spacing, we can 
neglect the inter-electrode block of the Hamiltonian and overlap matrix Hlr(rl) , Slr^rl^ (this also means that we 
neglect the direct tunneling between the two electrodes) and it is straightforward to solve G MM as: 





— {E + Smm 


— Hmm — T*l(E) — E^ 




EM 


= (E + Sml - 


- Hml)G°£l(E + Slm — 


Hlm), 


Eg(£) 


= (E+Smr - 


- Hmr)G r ^(E + Sum ~ 


Hrm) 


n 0;R 
^LL 


= {E+S LL - 


Hll)-\ 




n 0;R 
RR 


= (E + S RR - 


Hrr)' 1 . 





(3.11) 



Eq. ( |3.11 ) expresses the Green's function in the molecule in terms of the Hamiltonian matrix element in the same 
region, with the coupling to the left and right electrode included rigorously in terms of the self-energy operators Ef (E) 
and Eg (E). Note again that due to the short-range nature of the basis set, only the finite block of (*ll(rr) is needed 
for the calculation of T^, R dE) corresponding to the orbital basis in the left(right) electrode that have non-negligible 

overlap with the orbital basis in the extended molecule. So the calculation of G MM involves matrix operations only 
on finite matrices. 

The matrix self-energy operator can be taken as the matrix elements of a non-local operator in real space: 



= J Ja&mWW) ( 3 - 12 ) 

= y^(i^ + 5W_L;im — L\im) mn {E^~SlM ; nj — ^LM;7ij) 

n in 

= £[ / drcp*{r)[E + H]cp rn {r)]G£^ mn [ f ^ n (f)[E+ ~ B\^(f)\ 

mn J •> 
^ mn 

drd?4>*{r){[E + - H}G°£«(r,7)[E+ H\}<t>itf) 

with the equation for E^ taking the same form, .— . 

Using the analytic continuation rules of LangrethEJ, we can write down the corresponding "lesser" self-energy matrix 
as: 



= ~^^{E + SML-im — HML;im)G°£^. mn (E + ShM;nj — HhM-nj) (3.13) 

mn 

Since the left electrode in taken to be in thermodynamic equilibrium, we have G£'£{E) = i(G°£ R (E) — G°£^(E))f(E — 
and therefore: 

E< ;i .(£) = i(E« ;ij - E^.)/(E - it L ) = iT L . tij f{E - it L ) (3.14) 

where ^L;ij = ^L ij ~ ^L-ij anc ^ f(E — ^l) = 1+e {E-^ L )/kT is the Fermi distribution function. The equation for E^ 
follows by replacing L with R in the above equation. 



Similarly from the kinetic equation Eq. (2.8), we can write down the correlation function in the 
"extended molecule" in terms of the basis function <j>i as G < (r 1 r'\E) — Gfj(E)(f>i(f)(j)*(f l ) = 

/dr^^^G^)^)^^ After rearranging the order of summa- 

tion and integration and integrating over f" and f" , we get G< (r, r";E) = £ y {r)4>* (f) £ fe „ G ik ( E )^k n ( E ) G nj ( E ) ■ 
Comparing the coefficients of <pi(f)<j>*(f') 1 we obtain the following matrix equation: 

G<{E) = G R (E)E<(E)G A (E) (3.15) 
If there is no inelastic scattering due to the electron- vibronic coupling, we get: 
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£<(£) = E<(E) + £<(£) = iT L (E)f(E - fi L ) + iT R (E)f(E - m ) (3.16) 

since there is no "lesser" self-energy operator associated with V xc . And we can express the correlation function in 
terms of the distribution in each electrodes: 

G<(E) = i[G R {E)T l{E)G a (E)]f (E - Mi ) + i[G R {E)T R {E)G A {E)]f{E - m ) (3.17) 

where the products within the brackets are matrix products. Every physical observable of interest can be computed 
from the matrix correlation function Gfj. In particular, the current density is: 

J(r) = J dEJ{f-E) = 1/2 5Z / dEG ij( E ) i, im i v ' - V)&(f)$(f ) (3.18) 

The terminal current can be calculated numerically by integrating the current density J(r) over the boundary 
surface between the molecule and the electrodes or any cross sectional area in the "extended molecule" due to the 
current continuity in the "extended molecule" region. But often it is more useful to compute the terminal current 
directly from the matrix Green's function and the-jnatrix self-energy operators. This can be achieved by defining a 
current operator, as in Caroli et al. and Datta E2t 

J(f, r';E)= e/h[H{f)G < (r, f* '; E) - G< (7 , r; E)H (?)} (3.19) 

whose diagonal element gives the divergence of the current density: 

J(r, f; E) = V • J(r; E) (3.20) 

The current over a surface enclosing the "extended molecule" is written as: 

hot = J dE j dSV ■ J{r-E) 
= J dE J drI(r,f;E) 

= e/h JdE J df^[H(f)G^E!)M^^- G ^ E )M^U^ H ^ 
= e/h J dEj2[H Jt G<(E) - G<H 3l ] 

= e/h J dETr[HG<(E) - G<{E)H] (3.21) 

Now again we have transformed the integral over the coordinates to the matrix equation involving the Hamiltonian 
and correlation function matrices. From here on the usual derivation using the matrix notation, often used in the 
second-quantized form or with the semi-empirical Hamiltonian matrix, can be .carried through without change and 
we get the familiar final form for the terminal current in the matrix notation asEEj: 

h(R) = e/h J dETr{T L(R) [f(E - ii L{R) )A{E) + iG<{E)]} (3.22) 

where A(E) = i(G R (E) - G A (E)) and 

I = I J dETr[T L (E)G R (E)r R (E)G A (E)}[f(E - Ml ) - f(E - fx R )} (3.23) 

where 

r L (E)=i(X R (E)-[Z R (Etf), (3.24) 

T R (E) = i(E*(£0 - [££(£)]*) (3.25) 

Note that in order to represent the coupling to the electrodes as self-energy operators, it is essential that the "extended 
molecule" can be described by an effective single-particle theory. In other words, the Hamiltonian of the "extended 
molecule" in its second-quantized form can be diagonalized. If the Hamiltonian describing the the "extended molecule" 
contains product over four creation/annihilation operators, the simplified-description of the electrodes using self-energy 
operator breaks down and more complicated expressions are neededrTi 50 !. 

The above formulae are powerful formal results. Its practical application will depend on the existence of efficient 
algorithms for accurate evaluation of the-, Hamiltonian matrix elements over local orbital functions. In practice, a 
Gaussian-type orbital basis is often usecEa. 
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IV. SELF-CONSISTENT FORMULATION: DEVICE AT EQUILIBRIUM 



A. Model of the molecular device 



The present study considers a single molecule sandwiched between two semi-infinite metallic electrodes. In addition, 
there can be atomic-scale features on the metal surface. In the above formulation of our matrix Green's function 
approach, we have not specified the form of the Hamiltonian operator. The specific form of the Hamiltonian comes 
into play as we consider the self-consistent formulation. I— . .— . 

The analysis presented here is based on the spin-extension of Kohn-Sham density-functional (DF) theoryEZTE3 which 
requires the solution of effective single-particle, Schrodinger-like equations with spin-dependent potentials: 



H^1(r) = h~V 2 + V lon {r) + J ^ + V^rM^ = 

■-dp" 



00 - (4.D 



where a — a, (3 and 



p(r)=p a + p^ 

p"(r) = ^<K Q (f)| 2 , (4.2) 



Here nf is the occupation number of the spin orbital ip£ . For our system, the occupation of the eigenstate is governed 
by the metal Fermi-level Ep, so at low temperature, n\ — 0(ef — Ep ). In the valence-only calculations, the ionic 
potential Vi on {r) is represented by the non-local pseudopotential V ps {r 1 r l ), and the wavefuctions considered are the 
valence pseudo-orbitals. 

The spin-density functional theory (SDF) is the necessary generalization in the presence of magnetic field. If 
the external potential is only of electrostatic nature, as in our case due to the atomic nuclei, the DF formalism 
shows that it is possible to determine the system property using a functional that depends on the density alone 
and not the spin-densities. The main advantage of the SDF over the DF formalism is that the greater flexibility of 
the SDF formalism introduced by the spin dependence allows us to build more of the physics into the approximate 
functional. For example, the SDF formalism can give a reasonable description of the bond-breaking in molecules by 
allowing electrons of different spins to have different spatial density distribution. It also provides a better description 
of the open-shell molecules in satisfying the requirement (Hund's rule) that a state with a larger spin tends to be 



SDF formalism yields 
l-Hpur use of the SDF 



favored energetically. In addition, spin-orbit coupling effects can be included. As a result 
significantly better results for molecules and solids than the spin-unpolarized counterpart! 

formalism here follows that of the unrestricted Hartree-Fock method in molecular calculationsEJ in that we neglect 
the off-diagonal part of the spin density matrix, so the wavefunctions of the a- a»d /3- electrons are decoupled from 
each other except for the spin-dependence in the exchange-correlation potentialEZI(for molecules with even number 
of electrons and singlet spin states, the spin-unrestricted procedure usually leads to the same result as the spin- 
restricted one) . This breaks the rotational invariance in the spin space, and therefore cannot handle situations where 
the correlation between diffeteiit spin channels are important, e.g., the formation of a local magnetic moment at the 
molecule (the Kondo effect )IIjlZ3. These effects don't seem to be important, at least when the coupling between the 
molecule and the metal is strong. 

We will consider both the local spin-density approximation (LSDA) and its generalized-gradient approximation 
(GGA), the exchange-correlation energy of which takes the following form: 

E xc [n a ,np] = J d 3 rf(p a ,pj 3 ,j aa ,j a p,-f f3 j 3 ), 

laa = \Vn a \ 2 , (4.3) 
7a/3 = Vn a • Vnp, 
7,3,3 = IVn^l 2 
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B. The self-consistent formulation 



The starting point of computing the nonlinear transport characteristics of the molecular device, is then to com- 
pute the self-consistent charge and potential distribution of the device at equilibrium. This problem is equivalent 
to calculating the electronic structure of the metal-molecule-metal junction, which is the generalization of the famil- 
iar chemisorption problem where we consider the adsorption of a single isolated molecule onto one metal surface. 
Chemisorption is an important part of surface physics/chemistry a»d^y&s-Uptial theoretical studies exist r«|-tlic lit- 
erature, using the wasffifjipctiorjEa or Green's function. farmulatiorC^'OilZl'EaLa, and working in real spaceEJEZl or in 
finite basis expansionLfrO or combinations of the twcc3'L3L3. Our method is similar to these methods in many details. 
However, the uniqueness of our approach is that we use the finite basis expansion from the beginning, which allows us 
to separate the entire system naturally into different parts according to their geometrical arrangement, the treatment 
of which can then be improved independent of each other. This allows us to use existing, well-established techniques 
for treating the component systems — the molecules and the surface or the bulk — directly in the study of molecular 
devices with little change. 

As explained in the above, instead of solving the Kohn-Sham equation directly, we approximate the single-particle 
wavefunction by an expansion in a finite set of local orbital basis functions and solve the matrix Green's function 
Gfj (Z) as the solution of the matrix equation: 

J2(ZS ik -H° k )Gl j (Z) = S ij (4.4) 

k 

Here Z is an arbitrary complex energy. From the spectral representation relation: 

G^) = £f^, (4-5) 



we can obtain the density matrix , which is defined as 



^ = E^ c « ( 4 - 6 ) 



by performing the following contour integration in the complex energy plane: 

dZG^(Z), (4.7) 



1 

2ni 



c 



where the integr ation contour encloses the real energy axis from the energy of the lowest occupied state up to the 
Fermi energy i?^oE!l A typical integration contour is shown in Fig. (^). The advantage of integrating along the 
complex energy contour is that by moving away from the real energy axis, the sharp features in the density-of-states 
are smoothed, therefore allowing for accurate integration with a small integration meshEJ. 

The electronic charge density p a (r) , which is needed for the evaluation of H". , can be computed according to 

P%r) = Y J PlM^^ ( 4 - 8 ) 

We can also obtain the local density of states (DOS), defined as 

A* 

= E E c^bi^mE - ( 4 -9) 

from the following equation: 

n a (r 7 E) = -- lim V Imag[G?AE + i5)U i {r)(f) f Ar) (4.10) 



7T <5^0+ 



and the total density of states 
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n a (E)= ( d 3 rn a (r 7 E) = -- lim V Imag[GfJE + iS)]Sji = --Tr{Imag[GS]} (4.11) 

/ 7T S— >0+ * — ' 7T 

We can also obtain the projection of the total density of states into the molecular region as 

n a Mo i{E) = ~-Tr{Imag[GS]\ M oi} (4.12) 

where the trace is taken with respect to the orbital indices of the molecule. The transmission coefficient through the 
molecule is determined from: 

T{E) = Tr{T L G R T R G A }, 

r L =i(££-(E£)t), (4.13) 

r fl = i(sg-(sg)t). 

From here on, we avoid writing down explicitly the spin indices unless otherwise noted, with the understanding that 
we need to solve the Green's function corresponding to each spin direction and the total transmission is the summation 
over the two spin channels. 

The self-consistent computation is greatly simplified by realizing that due to the metallic screening within the 
electrode, the charge and potential perturbation induced by the adsorption of the molecule extends over only a fb 
region into the electrodes, practically only the region including the surface metal atoms closest to the molecuL 
The charge and potential distributions beyond this region are the same as that of the bimetallic interfaces without 
the molecules. It is this region — the molecule plus the perturbed surface atoms — that enters the above matrix 
Green's formulation and is called the "extended molecule" . The computation of the Green's function of the left and 
right electrodes then depends only on the charge and potential distribution of the bimetallic interfaces, which can 
be calculated separately and need to be calculated only once. Note that in this formulation we have neglected the 
long-range charge perturbation to the metal surface due to the presence of nonzero electric field in the "extended 
molecule" . This long-range charge perturbation could be important if we want to compute quantities that increase 
with distance such as the dipole moment associated with the adsorbed molecule. Since we are mostly interested in 
the electrostatic potential in the "extended molecule" region which is inversely,pcoportional to distance, no significant 
error will be introduced by neglecting these long-range charge perturbationsE3n3. 

The device characteristic is determined only by the electronic processes in the finite "extended molecule" region. 
In the MGF formalism, the effects of the contacts enter in two different ways: (1) as the infinite electron source and 
drain, the contacts inject electrons into and absorb electrons from the molecule, the occupation of the single-particle 
states within which is set by the electrochemical potential of the contacts. The self-energy operators describe this 
interaction between the molecular states with the continuum of states in the contacts, which of course depends on the 
band structure of the metal. (2) as the boundary condition, the charge and potential distribution in the "extended 
molecule" must join the charge and potential distribution deep within the interior of the contacts. Since the metallic 
electrodes are described well by the local-density-approximation, we assume only the exchange-correlation term in 
the "extended molecule" region may have a nonlocal form. So only the determination of the electrostatic potential is 
constrained by the long-range coulomb interaction such that it joins to the bulk value at regions beyond the "extended 
molecule" . The exchange-correlation potential depends only on the charge distribution in the "extended molecule" 
region, so it is irrelevant when we discuss the constraint of the boundary condition on potential imposed by the 
presence of the electrodes. 

Note that the unperturbed part of the contacts are the semi-infinite crystal with the perturbed surface atoms 

0' R 

removed. Its Green's function G£ L , RR -, can be calculated from that of the semi-infinite crystal using the "reduced 
space" idea of Williams, Feibelman and LangS, which is essentially a two-component version of our MGF equation 



(Eq. (3.10)). In this approach, we deal with the following two-by-two block matrix equation: 



E + S.LL{RR) ~ H LL ( RR } E + S L ( R ) P - i?L(_R)p \ I ^S;LL(RR) G S;L(R)P \ _ ( IS;LL{RR) 



E + S PL ( R ) - H PL ( R ) E + Spp - Hp P J I G ! f;PL(fl) Gg.pp J \ Is:PP 



(4.14) 



where we have separated the semi- infinite surface into two parts with P denoting the region of the semi-infinite surface 
that is included in the "extended molecule" and L(R) denoting the region of the semi-infinite surface corresponding 
to the unperturbed part of the left (right) electrode. From the above equation, we get: 

G LL(RR) = ( E+S LL(RR) ~ H LL ( RR ))~ l (4.15) 

_ riR. riR rr<R 

— U S;LL(RR) ~ u S;L(R)P\ U S;PP S U S;PL(R) 
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where again only matrix operations over finite matrices are involved due to the shore-range nature of the basis 
functions. 

The problem associated with the long-range electrostatic potential is commonly attacked using the scattering theory 
of electronic structure by writipa-dawn the Hamiltonian of the total system as a summation of an unperturbed part 
and a localized perturbationr 3 lr 9 H 76 H 7 % which is then solved using the Dyson equation, 

H = H° + AV eff = -iv 2 + V°(f) + AV eff , (4.16) 
G = G° + G°AV eff G (4.17) 
where the change in the effective potential is: 

AV eff [ P ] = AV lon + [ dV P ° (r>)] + V xc [ P (r)] - V xc [p°(f)] (4.18) 



In the above formulation, the "reference potential" V° and the "reference charge density" p° are those obtained from 
the self-consistent calculation of the bimetallic interfaces (without the molecule). In this way, AV e f* is nonzero only 
within the "extended molecule". AV' lon is the ionic potential due to the atomic ions in the isolated molecule plus 
any changes in the ionic potential that may be caused by the surface reconstruction and/or the modification of the 
molecular structure upon the adsorption onto the surfaces (any such structural relaxation is taken as occuring only 
in the "extended molecule" ) . 

Since the MGF formalism allows us to calculate the Green's function of the "extended molecule" given the potential 
within the same region, we want to transform the above equation into a form that better suits our purpose. Since 
the exchange-correlation potential in the "extended molecule" doesn't depend on the charge distribution outside this 
region, we don't need to include it in the choice of V°. This will also allow us to use different approximation schemes 
for the exchange-correlation potential in different part of the total system. Instead of calculating the change in the 
potential, we take a more constructive approach to get the potential directly and write down the Hamiltonian of the 
"extended molecule" as: 

H" = -iv 2 + V°(r) + V wn (f) - V ion '° + [ d 3 r + V° c (r), 



2 w W J 

h 

2 



_Iv 2 + V ext '°(r) + V lon {r) + f d 3 r'j^- + V° c (r), (4.19) 



where 

y«rt,0^ = v°(r) - V lonfi - I d 3 r . P °( r }. ( 4.20 ! 

J \r-r\ 

Here V ext '°(f) represents the electrostatic potential due to the charge (ionic and electronic) distribution in the contact 
region, .i.e., outside the "extended molecule". Since the contact region is charge neutral as a whole, V ext,0 (f) decays 
fast away from the contact region. If sufficient number of metal surface atoms are included into the "extended 
molecule", V ext '°(f) will be negligible within the molecule part. 

In principle, a separate self-consistent calculation of the bimetallic interface without the adsorption of the molecule 
needs to be performed^. For electrodes made from the same material, the results are practically the superposition 
of two separate surface calculations. For electrodes made from different materials, boundary conditions on the 
electrostatic potential will be imposed which line up the Fermi-levels of the two electrodesEil. Such a calculation 
will give us both the self-consistent charge and potential distribution in the region that we take to be beyond the 
range of admolecule perturbation, from which we can also get the corresponding surface Green's function G° LL and 
G RR . However, their exact value will depend on the particular model of the surface and the choice of the boundary 
between the perturbed and unperturbed surface region. Since the exact nature of the contact is almost never known 
in most experiments on molecular devices, and since the change in the electronic structure of the molecule is mainly 
determined by the local interaction between the molecule and the neighbor surface atoms that are included in the 
"extended molecule" , a good description of the contact can be obtained if we calculate the potential and the Green's 
function of the unperturbed part of the electrodes using values obtained from bulk calculations (for electrodes made 
from different materials, an additional linear potential term corresponding to the work function difference across the 
electrodes will need to be added to the bulk values for electrostatic potential calculations). This greatly simplifies 
the calculation while allowing meaningful quantitative comparison between theory and experiments. In particular, 
we note that satisfactory description of the charge and potential distribution in the bulk can be obtained without 
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performing fully seltconsistent calculations, for example, using the overlapping-atomic-potential mode£3, muffin- 
type approximations^ or tight-binding parameters. This approximate treatment of contacts can always be improved 
without affecting other parts of the c alculation. 



Given the Hamiltonian of Eq. ( 4.19 ), the matrix elements entering the Fock matrix in Eq. (3J3) can then be separated 
into the core, Coulomb and exchange-correlation matrix as usual in the molecular density-functional calculationsEil: 

H ij = H if re + Jl i + F xc;ij (4-21) 

where Hf° re is the summation of the kinetic-energy, the ionic potential and the "external potential" matrices, and 
Jij is the Coulomb matrix: 

kl 

PM = pti+pli ( 4 - 22 ) 

where we have used the conventiona l no tation for the electron repulsion integrals (ERI). For exchange-correlation 
energy of the functional form as Eq. ( |4.3| ), the exchange-correlation potential and their matrix elements are given by: 

W,/] = ^^, (4-23) 
Kc;ij = J V^*{r)H^r (4.24) 

and similarly for Vj? c and F^.^. 

At this point the outline of the calculation of the equilibrium property of molecular devices using the self-consistent 
matrix Green's function method is complete: 

(1) A judgment is made as to the set of atomic sites on the electrode surfaces which are included into the "extended 
molecule" . 

(2) A judgment is made as to the set of atomic sites on the electrode surfaces which are considered to be coupled to 
the "extended molecule" . 

(3) A basis set 4>i(r) is chosen for the atoms within the "extended molecule". 

(4) The Hamiltonian matrix of the unperturbed part of the electrodes and also the coupling between them and the 
"extended molecule" is determined from either the self-consistent bulk charge/potential distribution or from an ap- 
proximate construction such as the tight-binding model. The corresponding surface Green's function and self-energy 
matrix is then calculated at a preselected numerical energy integration mesh. 

(5) The "external potential" V ext '° and its matrix elements V^ xt '° are calculated either from a self-consistent sur- 
face/bulk calculation or by approximate construction. Also calculated are the matrix elements of the kinetic energy 
operator and the ionic potentials. Together they give the core Hamiltonian matrix. 

(6) An initial guess is made for the spin density matrix p a MM of the "extended molecule" . A natural choice for the 
initial guess is the density matrix of the free "extended molecule" . An even better one is the density matrix of the 
free "extended molecule" under the action of the "external potential" V ext '°. 

(7) Calculate the Coulomb and exchange-correlation part of the Fock matrix. For Gaussian-type orbital (GTO) basis, 
the Coulomb matrix elements can be calculated analytically. The existence of efficient algorithms for this operation 
is the main strength of GTO over other choices of local atomic orbital functions such as Slater- type orbitals (STO). 
The calculation of the exchange-correlation matrix elements must be performed numerically over the 3-d molecular 
volume d ue to the rather complicated form of the exchange-correlation potential (Eq. 



(8) Eq. (4.4) is then solved to obtain G a MM (Z) for each Z of the numerical integration mesh. The corresponding 



density- matrix p a MM is then recalculated via the numerical integration over the complex energy contour using Eq. 



(4.7) 



(9) Repeat step (7) and (8) until the input density or Fock matrix agrees with the output density or Fock matrix 
within a preset range. Due to the long-range nature of Coulomb interaction and the "heterogeneous" character of 
the "extended molecule" , strong oscillations often occur in the iteration process, therefor.e-|acceleration methods for 
self-consistent convergence are generally needed for the self-consistent process to convergeE3. 

(10) The electronic density of states, charge transfer, electrostatic potential and transmission coefficients are then 
calculated and analyzed. 

Note the above procedures are almost identical to that of calculating the electronic structure of the free "extended 
molecule" under some external potential V ext -°, the only exception being that at each iteration, instead of diagonalizing 
the Fock matrices, we integrate along the complex energy contour to obtain the density matrix for the next iteration 
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(see the program flow-chart of Fig. ||). In particular, given the density matrix pij, the computation of the Hamiltonian 
matrix is the same as that of the free "extended molecule" under the "external" potential V ext '°(f) which, like the 
self-energy operator, depends only on the charge and potential distributions obtained from the separate surface/bulk 
calculations and needs to be calculated only once. Since the above step is the computationally most demanding one in 
our self-consistent calculation, this allows us to take advantage of the full repertoire of molecular electronic structure 
calculations. In practice, this step can be replaced by call to any existing high-quality quantum-chemical software 
such as Gaussian 98E3. 



V. SELF-CONSISTENT FORMULATION: DEVICE OUT OF EQUILIBRIUM 

A. The meaning of the voltage drop 

Given the charge p° and the electrostatic potential V® s distribution of the device at equilibrium, applying an external 
bias will in general change the charge density and correspondingly the electrostatic potential (see Fig. ||), which are 
related through the Poisson equation: 

V 2 V es (f) = p(f), (5.1) 

V 2 K° s (r) = p°(r), (5.2) 

or 

V 2 5V es (r) =Sp(r) (5.3) 

where SV es — V es — V® g and dp = p — p° . The boundary condition on SV es is then: 

5V es {f)\ L = V L , 6V es (r)\ R = V R (5.4) 

where V^k) is the change of the electrostatic potential in the interior of the left (right) electrodes when current is 
flowing. V R — Vl — W could in general be different from V while p R — pl = eVE3. SV es describes the "voltage drop" 



across the molecule. The solution of Eq. (5.4) can be separated into two parts, which also makes its physical meaning 
clear: 



where Vbias(r) is the solution of 



5V es (f) = Vuas(r) + [ dr^^l (5.5) 
J \r-r\ 

V 2 V Mas (r)=0, 

Vbi aa (r)\ L = V L ,V Mas {r)\ R = V R (5.6) 

Note this is of the same form as we will obtain for the bare biased bimetallic interfaces. As is the case of the bare 
bimetallic interface, Vuas is linear. The final electrostatic potential is then given by: 

V es (r) = K° (r) + V blas (r) + [ tf^Qr (5.7) 

J \r-r\ 

where Sp is the charge redistribution in the "extended molecule" due to the applied bias. Written in this form, we see 
that Vbias doesn't depend on the property of the "extended molecule" , but only reflects the change in the boundary 
condition due to the applied bias and that is why we have used the subscript Vbias- The charge transfer 5p is the 
charge redistribution within the "extended molecule" induced by this "applied bias" Vbias- 

In this way, we can reach some qualitative conclusions about the voltage drop and the current- voltage characteristics 
of the molecular device. Obviously if the geometry of the molecule and the arrangement of the electrodes are such 
that the system is perfectly symmetric with respect to mirror reflection across the plane crossing the middle of the 
molecule, the current-voltage characteristic will be symmetric with respect to the bias. This can be seen from the 
above formula where changing the sign of V: V —> —V is equivalent to changing the coordinate z — > — z where we 
assume current flows along z axis. Since the choice of left or right is simply a matter of convention, the magnitude 
of the current flowing through the device can't change. Again in the symmetric case, since the charge neutrality is 
maintained in the "extended molecule" J dfSp(r) — 0, 5p(r) must be zero at the plane crossing the middle, so we 
expect the voltage drop will be close to the linear form in the middle of the molecule, and deviate from it as we move 
away from the middle towards the boundary set by the metal surface. 
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B. The self-consistent formulation 



Note the equilibrium electrostatic potential V° s is related to the equilibrium charge distribution p° through the 
following equation: 

V? s (r) = V extfi (r) + V lon {r) + / df^JQr (5.8) 

J \r-r\ 



From Eq. (5.7), the electrostatic potential at nonzero bias V e ° s is thus: 

V es (f) = V blas (f) + V ext <° (r) + V lon (r) + / (5.9) 

J \r-r\ 

Defining V ext (r) = Vbi as (r) + V ext: °(r), we ca n wr ite down the Hamiltonian of the molecular device at nonzero bias 
in the same form as that at equilibrium (Eq. (4.19)): 

H a = -iv 2 + V ext {f) + V lon {r) + f d 3 rj§^4-r + V x a c (r), (5.10) 



2 

Again V ext can be regarded as the electrostatic potential due to the charge distribution outside the "extended 
molecule" . The addition of Vuas reflects the change in the boundary condition in this electrostatic potentiaHO'H. 

For device out of equilibrium, the central quantity is the matrix correlation function Gf^ (E) , which is computed 
using the KKB equation: 

G<(E) = i[G R (E)T l{E)G a (E)]f (E - Mi ) + i[G R {E)T B {E)G A {E)]f{E - m ) (5.11) 

The self-consistent calculation then proceeds by computing the input density matrix to the next iteration from the 
correlation function computed in the current iteration: 

—G<(E)MrW 3 (r) (5-12) 



or in the matrix formulation: 



j^fm (5.13) 



The density matrix is nothing but the energy integ ration of the matrix correlation function. To see better the physical 
meaning of this, we divide both sides of Eq. ( 5.11 ) by: 

A(E) = i(G R (E) ~ G A (E)) = G R (E)(T L (E) + T R {E))G A (E) (5.14) 

and get: 

1l — G r (E)Tl(E)G a (E)/A, (5.15) 
7fl = G R (E)T R (E)G A (E)/A. 

where the division of matrices is defined such that = -g 12 -- Since the correlation function —iG < describes the 

number of electrons at energy E and the spectral function A describes the density of states at energy E, the above 
equation essentially says that the probability of the states at energy E being occupied in the molecule equals the 
probability of it being occupied in the left electrode multiplying the escape rate 7l from the left electrode to the 
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molecule plus the probability of it being occupied in the right electrode multiplying the escape rate 7^ from the right 
electrode to the molecule. 



Note the integration contour appearing in Eq. (5.13) is along the real energy axis rather than in the complex energy 
plane. Unlike the retarded Green's function, the correlation Green's function G < {E) is not analytic away from the 
real energy axis. This would have increased significantly the computational cost of the energy integration. Ho wever , 



for energy sufficiently lower than both /il and hr, say E < fJ-L(R) ~ lOfcT, we have f(E — (J-l(r)) ~ 1 an d Eq. ( 5.11 ) 
reduces to: 

G<{E) = i[G R (E)T L (E)G A (E)} + i[G R (E)T R {E)G A (E)] (5.16) 



Comparing with Eq. (5.14), we have: 

G< (E) = iA{E) = -2iImag[G R (E)} (5.17) 



Consequently the integration over the energy in Eq. ( 5.15 ) can be split into two parts: 



(Hi = i / dZGij(Z) + -L / dEGfj (E) (5.18) 

L ~ l JC L ~ % J E min 

where the first term represents integration along the same complex contour as that in Fig. ^ with the upper energy 
cutoff replaced by E m i n and the second term represents integration along real energy axis from E m i n to E max . Here 
E mm (max) is chosen such that for E < E min , f(E - (J, L {r)) w 1 and for E > E max J(E - (J. L ( R )) ~ 0. The integration 
along the real energy axis can be performed using fine integration grids since the range of integration is not much 
larger than eV. 

It is clear from the above discussion that the only difference between the self-consistent formul ation for device out 
of equilibrium and for device at equilibrium lies in the way of calculating the density matrix (Eq. ( |5.18| )). At nonzero 
bias, the retarded Green's function alone is no longer adequate for describing the observable property of the system. 
Instead, at each iteration, we need to calculate the correlation function from the retarded Green's function, the energy 
integration of which gives the input density matrix for the next iteration. All other procedu res re main the same. After 
self-consistency is achie ved, we can calculate the terminal current using Eq. (3.21) or Eq. ( |3.22| ) and also the current 



density using Eq. (3.18). As a result, calculating current transport in molecular electronic devices is no more difficult 



than solving a molecular chcmisorption problem, both conceptually and computationally. 



VI. CONCLUSION 



Modeling transport at the molecular scale requires a microscopic knowledge of the electronic structure of both 
the molecule and the surface which has to be considered in the context of an open system exchanging particles 
and energy with the external contacts. This complicates substantially the computational procedure and raises new 
questions about the theoretical basis. Although a first-principles theory of transport at the molecular scale doesn't 
yet exist, we have shown that the Non-equilibrium Green's Function Formalism of quantum transport combined with 
the density-functional theory of electronic structure provides a sound basis on which further works may be based. 

Since the electronic process in such molecular devices is mainly determined by the interaction between the molecule 
and the surface atoms closest to it, it is highly desirable to separate the treatment of the physical processes in this 
"extended molecule" region from that in the rest of the system since in general the exact atomic geometry of the contact 
is not known. In addition, as we have described, the treatment of the contact is best chosen to make the description 
of the physical processes in the "extended molecule" consistent. This is conveniently dealt with within the matrix 
Green's function approach using expansion in finite set of local orbital functions, where the effect of the contacts is 
taken into account using self-energy operators E^m. In addition, by introducing an "external potential" V ext which 
describes the electrostatic potential due to the charge distribution outside the "extended molecule" , we have shown 
rigorously that the electronic and transport property of the molecular device can be determined from the electronic 
processes in the "extended molecule" alone, given the knowledge of V ext and E^rm. This allows the highly accurate 
and efficient techniques developed for molecular electronic structure computation to be used for studying transport 
through molecules with little change. We believe such an approach will greatly accelerate theoretical/computational 
research on molecular electronic devices. 

Besides the computational advantage, the use of the local orbital basis sets and the techniques of molecular electronic 
structure theory will also greatly facilitate the interpretation of the result of computation using the language of 
qualitative molecular orbital theory which has provided the rationalization of the intuitive picture of bonding and 
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orbital interactions in moleculescJE3. Traditionally such discussion has relied on semi-empirical theories such as the 
Extended Hiickel type of theory. The advancement of density-functional theory in quantum chemistry has made 
the use of self-consistent Kohn-Sham orbitals in such molecular orbital theory highly desirableEirt3. For molecular 
electronic devices, the orbitals involved are those of the molecule and the surface atoms closest to it. Much of the 
physics can then be understood in terms of the orbital interactions after the effect of self- consistent charge transfer 
and potential distribution has been included. Such an approach has been taken recently by the present authors to 
study the. equilibrium property of the molecular device formed by the phenyldithiolate molecule bridging two gold 
contactsEj. Further work on nonlinear transport through molecules is under investigation^. 
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Notes added aftex^pzpof: After the completion of this manuscript, we became aware of two recent publications by 
Guo and coworkersEaK which use a formalism based on NEGF and density functional theory similar to that described 
here. However, there is significant difference in detail between our approach and theirs. 
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FIG. 1. Illustration of typical molecular devices. 
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FIG. 2. Integration contour in the complex energy plane. The cutoff energy Eo should be below the lowest occupied states 
of the system considered. 
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FIG. 3. Program flow chart of the self-consistent Matrix Green's function (MGF) approach to molecular device calculation. 
Note that besides the introduction of V ex t and the self-energy operator ^l(r) > the only difference with the conventional molecular 
electronic structure calculation lies in the calculation of the density-matrix given the molecular Fock matrix. See the discussion 
in the main text. 
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FIG. 4. (a) Schematic illustration of the electrostatic potential profile V® 3 in the absence of applied bias. V" e ° s joins to the 
bulk values in the interior of the contacts. The exact value of V^ s should be obtained from the self-consistent calculation of the 
metal-molecule-metal junction as described in the main text, (b) Schematic illustration of the electrostatic potential profile 
V ss in the presence of applied bias V. Note the changes in the boundary condition of V ea under bias, (c) Difference in the 
electrostatic potential 8V es = V es — V° 3 in the presence and absence of applied bias V. Full line shows the linear "applied bias" 
Vbias wh ile dotted line shows the true SV es including the effect of the self-consistent charge redistribution Sp under bias (see 



Eq. (5.5) and discussions thererein). 
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